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We have used a density-functional-based tight-binding method in order to create structural models 
of the canonical chalcogenide glass, amorphous (a-)As2S3. The models range from one containing 
defects that are both chemical (homopolar bonds) and topological (valence-alternation pairs) in 
nature to one that is defect-free (stoichiometric). The structural, vibrational and electronic prop- 
erties of the simulated models are in good agreement with experimental data where available. The 
electronic densities of states obtained for all models show clean optical band gaps. A certain degree 
of electron-state localization at the band edges is observed for all models, which suggests that pho- 
toinduced phenomena in chalcogenide glasses may not necessarily be attributed to the excitation of 
defects of only one particular kind. 



I. INTRODUCTION 

Amorphous chalcogenides (particularly the sulfides, se- 
lenides and tellurides) exhibit intriguing physical prop- 
erties that are not observed in their crystalline counter- 
parts. Some of these unusual properties are extensively 
used in electronic and photonic devicesi and there are 
many potential applications. 

Perhaps the most interesting opto-electronic behavior 
of these materials is the metastable structural changes 
resulting from the absorption of near-bandgap lights. 
The microscopic changes in the atomic structure involved 
are not generally observable directly, but they are re- 
flected in measurable optical, electronic, and mechani- 
cal properties^. A fully consistent microscopic theory 
accounting for the opto-electronic behavior in chalco- 
genide glasses, however, is still lacking. It is therefore 
of great interest to employ computer simulations in or- 
der to generate and study structural models of such ma- 
terials. Provided that these computer-generated mod- 
els compare well with available experimental data, they 
could allow one to monitor the photo-induced structural 
changes in the greatest possible detail at the microscopic 
level. Needless to say, such detailed information is not 
generally presently available from experiments. 

In order to study opto-electronic effects, one needs to 
perform quantum-mechanical calculations that are very 
time consuming. This imposes a severe constraint on the 
accessible system size of the simulated models. It is pos- 
sible, however, that small samples may be sufficient to 
capture much of the interesting photoinduced behavior 
due to a high degree of localization of the photo-excited 
electron-hole pairs^. A few ab initio studies of amor- 
phous chalcogenides, paying particular attention to elec- 
tronic properties, have been performed by Drabold and 
co-workers (see, e.g., Refs.QHI^)- Although their results 
agree well with experimental data, the defect concentra- 



tion in these structural models is much greater than is 
estimated in experiments. Some of the possible reasons 
for this are the level of approximation in the ab initio 
approach employed in these studies and/or the model 
preparation history, which may have resulted in too many 
quenched-in defects due to very rapid cooling of the sam- 
ple from the liquid state. 



In this paper, we use a density-functional-based tight- 
binding (DFTB) method^*^ in order to generate and ana- 
lyze several models of amorphous diarsenic trisulphide 
(a-As2S3) with a controlled and systematic change in 
the defect concentration. Apart from being widely re- 
garded as the canonical chalcogenide glass, this partic- 
ular material was chosen for analysis for the following 
reasons. First, we have obtained reliable high-quality 
neutron-scattering structural data for this materialiii^. 
Second, the properties of a-As2S3 are expected to be sim- 
ilar to those of a-As2Se3, amorphous diarsenic triselenide, 
that have been studied theoretically in Ref. d, and it is 
of interest to corroborate this similarity in independent 
simulations. Third, the necessary input DFTB data for 
sulphur have been previously generated and extensively 
testedi^. 



Although the DFTB method is semiempirical, it al- 
lows one to improve upon the standard tight-binding de- 
scription of interatomic interactions by including a DFT- 
based self-consistent second order in charge fluctuation 
(SCC) correction to the total energy. The flexibility in 
choosing the desired accuracy while computing the in- 
teratomic forces brings about the possibility to perform 
much faster calculations when high precision is not re- 
quired, and refine the result if needed. 
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II. METHODOLOGY 



A. DFTB 



The SCC-DFTB model is derived from density- 
functional theory (DFT) by a second-order expansion 
of the DFT total energy functional with respect to the 
charge-density fluctuations 5n' = 5n{r') around a given 
reference density tiq = no(r')- 
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where / df' is expressed by /'. Here, = H[no] is the 
effective Kohn-Sham Hamiltonian evaluated at the ref- 
erence density and the are Kohn-Sham orbitals. E^c 
and Vxc are the exchange-correlation energy and poten- 
tial, respectively and En is the core-core repulsion energy. 

To derive the total energy of the SCC-DFTB method, 
the energy contributions in Eq. JQ) are further subjected 
to the following approximations: 

1) The Hamiltonian matrix elements {'^i\H'^\'^i) are rep- 
resented in a minimal basis of confined, pseudoatomic 
orbitals 6,, , 



To determine the basis functions 0^ , we solve the atomic 
DFT problem by adding an additional harmonic poten- 
tial (7^)^ to confine the basis functions^^. The Hamilto- 
nian matrix elements in this LCAO basis, H^^, are then 
calculated as follows. The diagonal elements i/f?,, are 
taken to be the atomic eigenvalues and the non-diagonal 
elements H^^ are calculated in a two-center approxima- 
tion: 

which are tabulated, together with the overlap matrix 
elements S^i, with respect to the interatomic distance 
Raf3- Veff is the is the effective Kohn-Sham potential 
and n° are the densities of the neutral atoms a. 
2) The charge-density fluctuations Sn are written as a 
superposition of atomic contributions dua, 



6n — 6r 



which are approximated by the charge fluctuations at the 
atoms a, Aqa ^ ^ Qa- 9a is number of electrons 
of the neutral atom a and the qa are determined from 



a MuUiken-charge analysis. The second derivative of the 
total energy in Eq. ^ is approximated by a function 
Jaf3, whose functional form for a =/= (3 is determined an- 
alytically from the Coulomb-interaction of two spherical 
charge distributions, located at Ra and Rp. For a = /3 it 
represents the electron-electron self-interaction on atom 
a. 

3) The remaining terms in Eq. (Q, En and the energy 
contributions, which depend on uq only, are collected in 
a single energy contribution E^^^p. Erep is then approxi- 
mated as a sum of short-range repulsive potentials, 

Erep — ^ U[Ral3], 

which depend on the interatomic distances Ra/3- 

With these definitions and approximations, the SCC- 
DFTB total energy finally reads: 



(2) 



a/3 



Applying the variational principle to the energy func- 
tional ||2J), one obtains the corresponding Kohn-Sham 
equations: 

^Ci.j(i/^i. - e.j5^i.) = 0, y i (3) 

u 

^ C 

which have to be solved iteratively for the wavefunction 
expansion coefficients c^, since the Hamiltonian matrix 
elements depend on the due to the MuUiken charges. 
Analytic first derivatives for the calculation of inter- 
atomic forces are readily obtained, and second deriva- 
tives of the energy with respect to atomic positions are 
calculated numerically. 

The repulsive pair potentials U[Raf3] are constructed 
by subtracting the DFT total energy from the SCC- 
DFTB electronic energy (first two terms on the right- 
hand side of eq.lO) with respect to the bond distance 
Rai3 for a small set of suitable reference systems. 

To sum up, in order to parameterize the method for a 
new element, the following steps have to be taken. First, 
DFT calculations have to be performed for the neutral 
atom to determine the LCAO basis functions 0^ and the 
reference densities nj^. Here the confinement radius can 
in principle be chosen different for the density (rp) and 
each type of atomic orbital {rQ^''^). We usually take tq 
to be the same for s- and p- functions. In a minimal basis, 
this yields a total number of two adjustable parameters 
for elements in the first and second rows, while there 
are three if d-functions are included. After this, the dif- 
ferent matrix elements can be calculated and the pair 
potentials UlRap] are obtained as stated above for every 
combination of the new element with the ones already 
parameterized. 
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In this study, we used the same tabulated data set for 
sulphur-sulphur interactions as in Ref. m and the As- 
As, As-S and S-As data sets were generated according 
to the procedure outlined above. The confinement radii 
for the As pseudoatomic densities (r^ = 9.8 a.u. and 
= 4.5 a.u.), as well as the As- As repulsive pair-potential 
were determined in accord with other ongoing efforts re- 
lated to GaAs systems. The cage-like AS4S6 molecule 
was used to calculate the As-S repulsive pair potential, 
so that after finding the minimum energy configuration of 
the molecule in all-electron DFT-LDA calculations using 
the NRLMOL programi^, a regular scaling of the As- 
S bond-lengths was performed, keeping the overall T^ 
symmetry. Then the acquired potentials were tested on 
other clusters, such as AS2S and AS4S4 molecules, with 
an overall good agreement of the binding energies and 
configurations between SCC-DFTB and the reference all- 
electron DFT-LDA NRLMOL results. When these data 
sets were used in order to optimize the geometry of the 
crystal structure of orpiment {C-AS2S3) in SCC-DFTB, 
the agreement with the experimental structur*^^ was 
within 2%. 



B. Preparation of structural models 

In experiments, bulk glasses are usually prepared from 
the melt by rapidly cooling (quenching) the sample. Al- 
though it appears impossible to achieve experimentally 
realistic cooling rates in molecular-dynamics computer 
simulations, some empirical procedures result in mod- 
els that can be in good agreement with experiments. In 
order to prepare realistic models of a-AsoSa, we use an 
algorithm akin to that used, e.g., in Refs. ItIIMITI 

The structural model of a-As2S3 was obtained in the 
course of an NVT (constant number of particles, N , vol- 
ume, V , and temperature, T) molecular-dynamics simu- 
lation with periodic boundary conditions. Since we are 
not interested in statistical properties of thermal fluctu- 
ations, the temperature was controlled simply by scaling 
the velocities of the constituent particles every few time 
steps with the time intervals between the scalings taken 
randomly with a mean value of 10 time steps. We used a 
time step of 100 a.u. « 2.4 fs (1 a.u. = 2.4189x10"" s) 
and the Verlet algorithm in order to integrate the equa- 
tions of motion. 

The starting configuration was a system of 200 (80 ar- 
senic and 120 sulphur) atoms in a cubic supercell, with a 
side length of 17.25 A, obtained by rescaling a crystalline 
configuration of orpiment (monoclinic, space group 14, 
P12i/nl) with 1 x2x5 20-atom unit cells. The crystalline 
coordinates were obtained from the Inorganic Crystal 
Structure Database and correspond to those reported in 
Ref. In order to obtain a cubic supercell, we approx- 
imated the monoclinic unit cell of orpiment by an or- 
thorhombic one simply by neglecting the small deviation 
of the angle (3 = 90.68° from the right angle. Then we 
used the experimental glass density of p = 3.186 g/cvp? 



from Ref. in order to obtain the side length of the 
cubic supercell L — {N/ pY^^ =17.25 A and the coordi- 
nates of atoms from the cuboid with dimensions 11.48 
by 19.15 by 21.28 A were scaled by 1.5, 0.9, and 0.81 
in the x, y, and z directions, respectively. The use of 
the crystalline initial configuration gives the correct sto- 
ichiometry and the rescaling of a non-cubic supercell in 
order to obtain a cubic one serves the goal of eliminating 
possible anisotropies in physical properties. Note that in 
Ref.lll, the authors cut a cubic supercell from a crystalline 
phase of As2Se3 isostructural to orpiment, thus achieving 
the above two goals, but disrupting the periodicity with 
respect to the periodic boundary conditions. After thor- 
ough equilibration in the simulated liquid state, however, 
the use of either of these two prescriptions should lead to 
models with similar statistical properties. 

The initial configuration was melted and equilibrated 
first at T = 3000 K for 3 ps, and then the resulting con- 
figuration was allowed to equilibrate at T = 1000 K for 12 
ps. The equilibration criterion used was the convergence 
of the average potential energy to a constant value. At 
such high temperatures, the accuracy of the calculations 
appears to be least significant for the subsequent genera- 
tion of low-temperature structural models. Therefore, at 
this stage, a minimal basis set of only s and p orbitals on 
both the As and S atoms was used, and the tight-binding 
scheme of Sec. Ill Al was used without the self-consistent 
charge (SCC) correction in order to speed up the cal- 
culations. It was necessary, however, to use the SCC 
correction during the initial 1.5 ps of the T = 3000 K 
run, since the large forces resulting from the presence 
of small interatomic distances in the distorted starting 
configuration otherwise led to numerical instability. 

While preparing computer-generated models of glasses, 
it is customary to mimic real experiments by reducing the 
temperature over time intervals whose length, however, is 
limited by the available computer time. We found that, 
due to the unrealistically small length of such time in- 
tervals, this approach is rather impractical. Instead, we 
used the available computer time to perform an anneal- 
ing run at one fixed temperature which is low enough 
for the process of the bond-network formation to be ac- 
tivated and high enough for the topologically connected 
network to grow sufficiently rapidly. First, we performed 
a run corresponding to 6 ps at T = 700 K with the 
SCC correction and the minimal (sp) basis. Keeping in 
mind that the simulation was done at constant volume, 
this temperature was chosen to be somewhat above the 
melting temperature of orpiment at atmospheric pres- 
sure {Tra — 592 K according to Ref. IT^ . We empirically 
found that annealing the configuration for the following 
6 ps at a higher temperature of T = 800 K slightly im- 
proved the quality of the network by increasing the frac- 
tion of heteropolar bonds in the model. During these 
two runs, we used a smaller time step of 50 a.u. (1.2 
fs). Finally, the temperature was nearly instantaneously 
reduced to T = 300 K, quenching the system within a 
metastable basin on the potential-energy hypersurface. 
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The resulting model (model 1 in the following) remained 
stable while, during a run corresponding to 120 ps, 500 
configurations separated by 10 time steps of 100 a.u. (24 
fs) were stored for subsequent analysis. In the last run, 
we used the SCC correction and increased the basis set 
by including the d orbitals for sulphur atoms. This ba- 
sis set extension provided a major impovement in the 
description of hypervalent sulphur molecules ^'^ as well as 
silicon-oxygen compoundsSS. While d orbitals on sulphur 
atoms give noticeable improvement, it appears that they 
are less important for arsenic and we restrict ourselves to 
including only s and p orbitals for the As atoms in order 
to speed up the calculations. 

Model 1 contains three topologically identical coordi- 
nation defects, namely intimate valence alternation pairs 
(IVAPs), where a singly coordinated sulphur atom is at- 
tached to a three-fold coordinated arsenic atom, thus in- 
creasing the coordination number of this arsenic atom to 
four. Apart from the IVAPs, the amorphous network is 
topologically ideal, in the sense that each sulphur atom is 
bonded to two neighbors and each arsenic atom is bonded 
to three. There is, however, a certain degree of chemi- 
cal disorder in this system which manifests itself in the 
presence of nine As- As and six S-S homopolar bonds. 

In the context of photoinduced metastability, a great 
deal of significance is attributed to the presence of topo- 
logical and/or chemical defects^. It is therefore impera- 
tive to create models both with and without such defects 
in a theoretical investigation that attempts to be con- 
clusive. We produced additional models by "surgically" 
removing the IVAP defects and homopolar bonds from 
model 1. Model 2, which does not contain any topologi- 
cal defects, was obtained by removing the three singly co- 
ordinated sulphur atoms from model 1 and rescaling this 
197-atom model to the original density. This procedure 
did not affect the stability of the amorphous network. 

In order to eliminate the chemical defects, we itera- 
tively applied the following algorithm that utilizes the 
ideasSi used to create models of binary amorphous solids 
(e.g. a-Si02) from one-component continuous random 
networks (e.g. a-Si)^^. First, a sulphur atom was in- 
serted in the middle of each As- As homopolar bond. Sec- 
ond, each S-S bond was replaced by a single sulphur atom 
located at its mid-point so that each local As-S-S-As con- 
figuration turned into As-S-As. Third, the distance be- 
tween each newly introduced S atom and its two nearest 
arsenic atoms in the newly created As-S-As units was re- 
duced in order to increase the bonding character of the 
As-S bonds stretched by the above manipulation. We set 
a constraint on the length of the modified As-S bonds so 
that it did not exceed 2.5 A. Fourth, the modified con- 
figuration was relaxed in an MD run a,t T — 300 K for 
2.5 ps, until the potential energy reached a plateau. Af- 
ter the first iteration, the 200-atom sample (that we call 
model 3 in the following) contained only one As- As and 
one S-S bond that were spatially well separated (the min- 
imum As-S distance among these four atoms was 4.6 A). 
Only two iterations were sufficient in order to obtain a 



model with all-heteropolar bonds, which we refer to as 
model 4 in the following. The defect statistics for models 
1-4 are summarized in Tabled 



C. Data analysis 

A common way to assess the quality of a structural 
model is to compare experimental and calculated static 
structure factors. We have calculated the structure factor 
S{Q) by Fourier transforming the radial pair-correlation 
function g{r) (also often called the pair- or radial distri- 
bution function) defined as (see, e.g., Refs. I2^l2^l25l) : 



V 



47rr2iV2 




(4) 



where the sum is over all pairs of atoms in the sample 
of volume V separated by distance Vij, N is the total 
number of atoms and the angular brackets denote an en- 
semble average. In the case of neutron scattering, that 
is of interest here, bi is the coherent scattering length 
of atom i and (6) is the average scattering length. This 
function gives the probability of finding a pair of atoms a 
distance r apart, relative to the probability expected for 
a completely random distribution of atoms at the same 
density. For a binary alloy, e.g. AS2S3, it is of interest to 
decompose g{r) in terms of the partial pair-correlation 
functions gap{r): 



9 



(5) 



where the double sum is over atomic types and ba = 
Cabal {b), with Ca = Na/N being the atomic fraction of 
a atoms. From Eq. Q), it follows that 



V 



\ia^]f3 



(6) 



where the index ia runs over a-type atoms only. The 
values of the scattering lengths used here were 6as = 
6.58 fm and bs = 2.847 fm (see, e.g., Ref.l2^. In practice, 
we use a standard algorithm, where the (^-function in 
Eq. © is replaced by a function which is non-zero in a 
small range of separations, and a histogram is compiled 
of all pair separations falling within each such range (see 
e.g. Ref. I27I). Analogously, the bond-angle distribution 
function can be calculated as a histogram of all bond 
angles in the system. 

While bond-angle distributions provide information on 
the short-range order of an amorphous material, ring 
statistics have been generally used as a measure of the 
medium-range order. An n-membered ring is a closed 
loop with n atoms (or bonds). Here, we count only the 
shortest-path (irreducible) ringsSS, i.e. those which do 
not have "shortcuts" across them. In order to identify 
such rings we use the algorithm due to Franzblau^^, as 
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implemented in the program "statix" by Jungnickel22i^. 
The basic idea is to travel along the network paths (bond 
chains) containing a tagged atom and identify all the 
rings of length up to a given maximum. For each ring, 
it can then be verified whether it is an irreducible one 
(containing no cross links). 

Experiments also provide information on the vibra- 
tional and electronic densities of states (VDOS and 
EDOS) which can be compared with the results obtained 
from simulation. In addition to this, simulations allow 
one to assess the degree of localization of the vibrational 
and electronic eigenstates. We compute both the vibra- 
tional and electronic densities of states by using the fol- 
lowing definition: 



g{uj) = C^S{uj -uj„), 



(7) 



where the constant C is determined by normalization, aj„ 
are eigenfrcquencies of the Hessian (dynamical) matrix 
of an energy-minimum configuration in the case of the 
VDOS or Kohn-Sham eigenfrcquencies (or energies) cor- 
responding to this configuration in the case of the EDOS, 
and the sum is over all eigenstates. In practice, in order 
to obtain a smooth representaion of g{Lu), the delta func- 
tion in Eq. (jT)) is replaced by a Gaussian function centered 

at UJn- 

In contrast to the eigenstates of a perfect crystal, that 
extend over the entire sample, some eigenstates in dis- 
ordered solids are localized at relatively small groups 
of atoms. The degree of localization can be quanti- 
fied by the inverse participation ratio that is defined in 
terms of dynamical-matrix eigenmodes or MuUiken par- 
tial charges, for vibrational or electronic excitations re- 
spectively. For a vibrational mode n, the inverse partic- 
ipation ratio can be defined^ as 



Model 



No. of atoms 
No. of As- As bonds 
No. of S-S bonds 
No. of IVAPs 
Total No. of defects 



200 197 200 200 
9 9 10 
6 6 10 
3 
18 15 2 



TABLE I: Defect statistics for models 1-4. 

where the contribution to state n from atomic site i, 
can be expressed in terms of the wavefunction coefficients 
in the tight-binding basis c^„ and the elements of the 
overlap matrix S : 



(10) 



Here, in the double sum, the index /i runs only over the 
atomic orbitals located on atom i and the index v goes 
over all orbitals. As in the case of the vibrational inverse 
participation ratio, the electronic is equal to 1/N for 
a totally delocalized mode and approaches unity with in- 
creasing degree of localization. The partial charges g-""* 
allow one to detect on which atoms most of the charge is 
localized for a particular eigenstate. By refraining from 
summing over the atomic orbitals, i.e. over /i in Eq. I|1U|I . 
one can identify the type of the atomic orbitals, e.g. s 
or p, most actively participating in an eigenstate. Anal- 
ogously, the local EDOS for a particular orbital type can 
be obtained via the following expression : 

g^{uj) ^ C'^6{lO - UJn)'^S^uCfj,nCvn (H) 
n v 




(8) 



When the displacement eigenvectors ' , n 

1,2,...,3A^, are normalized to unity (X^il'^i"''^ ~ 1)' 
= 1 for a mode totally localized at one atom and 
— 1/N for a completely extended mode, such as a 
rigid-body displacement. 

In the case of the electronic properties, the linear com- 
bination of atomic orbitals (LCAO) concept employed in 
the DFTB program allows one to separate the contribu- 
tions from individual atomic sites and orbitals to the total 
charge for a particular eigenstate, and to decompose the 
total EDOS in terms of the local electronic densities of 
states (LEDOS). Using the MuUiken population analysis 
(see, e.g., Ref. 1321), the inverse participation ratio for an 
electronic state n can be written as^ 



Pn = 



N 

E 

i=l 



,(")|2 



(9) 



III. RESULTS 



A. Structure 



By using the construction method described in 
Sec. Ill Bl we obtained four models of a-As2S3 which 
are distinguished by the presence and concentration of 
topological and chemical defects (see Table IJ). We also 
created 60- and 100-atom models with similar concen- 
trations of homopolar bonds as in models 1 and 2, and 
without coordination defects. 

Fig. nja) shows that the pair-correlation function 
(PCF) corresponding to model 1 compares well with 
two independent neutron-scattering experimental results. 
The discrepancy between the two experimental PCFs al- 
lows one to estimate the uncertainty in the experimental 
data. The main difference between the experimental and 
simulated results is in the height of the first peak. Since 
the experimental PCFs are obtained by Fourier trans- 
forming the measured static structure factor, where the 
large-Q oscillations are damped by applying a window 
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FIG. 1: (a) Total pair-correlation functions for model 1 
and the neutron-diffraction experiments 1 (Ref. 12) and 2 
(Ref. Il^ . (b) Total pair-correlation functions for models 1 
(the same as in (a)) and 4. (c) Partial pair-correlation func- 
tions for model 1. gs-sir) and gAs~s(r) are shifted upwards 
by 3 and 6 units, respectively. 



FIG. 2: (a) Reduced structure factors (interference func- 
tions ] fo r model 1 and the neutron-diffraction experiments 1 
('Ref. ll2t) and 2 fRef. llSl) . (b) Total interference functions for 
models 1 (the same as in (a)) and 4. (c) Partial interference 
functions for models 1 and 4. The functions corresponding 
to S-S and As-S correlations are shifted upwards by 3 and 6 
units, respectively. 



function, this reduces the height of the first peak in g{r) 
and also broadens its width. 

Although the PCFs corresponding to models 2-4 are 
quite similar to that for model 1 (which is why we do 
not show here the PCFs for models 2 and 3), there is one 
conspicuous distinction in the shape of the first peak that 
depends on whether or not the system contains homopo- 
lar bonds (see Fig.^b)). While this peak is symmetric 
for model 4 with all heteropolar bonds, there is a shoulder 
on either side of the peak in the PCF for model 1. From 



Fig. njc) , it is seen that the shoulders in the first peak 
of the total PCF originate from the homopolar As-As 
and S-S bonds which produce small peaks in the respec- 
tive partial PCFs at this position. It is remarkable that 
the calculated PCF for model 1 virtually reproduces the 
right-hand side of the first peak in the PCF from exper- 
iment 1 (see Fig. Wi''^)) ■ This result supports one of the 
conclusions of Ref. |13 (experiment 1) that the low- and 
high-r sides of the base of the first peak in the PCF can 
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FIG. 3: Close-up of the parts of F{Q) — Q{S{Q) — 1) most sensitive to changes in the local structure. The FSDP region for 
(a) the 60- and 100-atom models, model 1 and experiment 1, and (b) models 1-4, respectively. The region near Q ^ 7.5 A"^ 
for (c) the 60- and 100-atom models, model 1 and experiment 1, and (d) models 1-4, respectively. 



Sample 


rAs-s [A] 


^As-S 


^S-As 


experiment 1 


2.27 


2.8 


1.8 


model 1 


2.279 


2.81 


1.88 



TABLE 11: Values of average bond length and coordination 
numbers for the first coordination shell found from peak fit- 
ting (experiment 1) and direct calculation (model 1). 



be ascribed to S-S and As-As bonds, respectively. Ta- 
ble ^further demonstrates that the agreement between 
the structural characteristics of model 1 and experiment 
1 is very good and, in particular, that the As-S coordi- 
nation numbers for these two samples are practically the 
same. The larger discrepancy in the S-As coordination 
number and the low-r side of the first peak in the PCF 
arises from the relatively small system size and the fact 
that the numbers of As-As (nine) and S-S (six) homopo- 
lar bonds are not equal to each other in model 1. The 
position of the first peak in the partial PCF gs~s{r) also 
plays a role here. As is seen from Fig. ^c), this peak 
shifts towards the low-r end when d orbitals on sulphur 
atoms are included into the basis set. When a.T ~ 300 K 
run is performed without the d orbitals in the basis set, 
the position of the first peak in gs~s coincides with that 



of the first peak in ^as-Sj a-nd the shoulder on the low-?' 
side of the simulated total g{r) is not seen at all. 

It is instructive to compare the structural data also in 
Q-space, as this often emphasizes features that are not 
obvious in an r-space representation. Fig. |2Ia) shows 
the reduced structure factors (or interference functions) 
F{Q) = QiS{Q) - 1), related to the PCFs in Fig. [Ja) 
by a Fourier transform. Again, the agreement between 
model 1 and experiment 1 is very good. From the rate 
of decay of F{Q)iroin experiment 2, it is apparent that 
the data in Ref. [13 are reported for F{Q) multiplied by 
a window function. 

F{Q) for model 1 exhibits a first-sharp diffraction peak 
(FSDP) at about the same position, Q « 1.5 A~^, as 
found in the experiments. The magnitude of this peak 
depends on the system size and is expected to increase 
for a larger model. This statement is supported by the 
observation of this tendency as the peak develops in our 
60-, 100-, and 200-atom models (see Fig.E^a)). Fig.OIb) 
demonstrates a systematic displacement in the position 
of the FSDP towards the high-Q end as the number of 
defects is reduced from model 1 to model 4. In Ref. 
it was observed that the FSDP in X-ray diffraction in- 
tensity curves for an a-As2S3 thick (8 fim) film reversibly 
moved towards the low- and high-Q end upon illumina- 
tion and annealing, respectively. The result in Fig. Olb) 
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FIG. 4; Bond-angle distributions. Solid lines in (a)-(e) and dashed lines in (c) and (d) correspond to models 1 and 4 respectively. 
The height of vertical lines topped by circles is proportional to the number of distinct angles within one degree in the crystal 
structure of orpiment. Solid vertical lines and circles correspond to the experimental data from the crystallographic database 
and Ref. IT^ and the dashed vertical lines and open circles correspond to the structure optimized by the DFTB method. 



is consistent with the above observation, if we suppose 
that the defect concentration is reversibly increased and 
decreased upon illumination and annealing, respectively. 
The height of the measured^ peak, however, increased 
as its position wavenumber, Qi, decreased after illumina- 
tion, while in the simulations we observed the opposite 
tendency for the height of the FSDP to increase upon 
elimination of defects and increasing Qi. A possible rea- 
son for the different behavior of the height of the FSDP is 
that the simulations of the bulk a-As2S3 were performed 
at constant volume, while the experiment was done for 
an amorphous film at atmospheric pressure. 

Another distinction between the different interfer- 
ence functions presented in Fig. |21 is seen in the range 
7 ^0^8 A~^. F{Q) appears to be very sen- 
sitive to structural differences in this particular range of 
wavenumbers, as is apparent from Fig.|2Ib), where the in- 
terference function for model 1 is compared with that for 
model 4, and from Fig. |2Ic,d), where F{Q) is magnified 
in this Q-interval, and the differences between the curves 
are prominent. Although the differences between the par- 
tial interference functions for models 1 and 4 (shown in 
Fig.|3c)) are each rather subtle in this Q-interval, they 
become more pronounced when combined into the total 
F{Q) (see Fig. I^b)). The small peak, which is seen in 
F{Q) for model 1 at Q ~ 7.5 A~^, diminishes from model 
1 to 4, so that it is practically not seen in the case of the 
stoichiometric model 4 (see Fig.l^Jd)). The presence and 
position of this small peak may be attributed to the pres- 
ence and spatial distribution, respectively, of homopolar 
bonds in the system. Similar peaks are seen in F{Q) at 
about Q = 7.5 A"-'^ for both experiments mentioned here 
(Fig. [^LlL R-efs. ^2 18) and in the experiment reported 
in Ref. fill While, in all these independent experiments, 
these peaks virtually coincide, in the different models pre- 
sented here they do not agree so well (see Fig.lSJc)), as, 
perhaps, can be expected in the case of small system 
sizes. 



Another structural characteristic that is of interest, 
and that is easily accessible in computer simulations, 
is the bond-angle distribution (see Fig. Although 
we have no experimental data for this distribution, the 
known structure of a corresponding crystal can serve as 
a guide in assessing the quality of our models — many 
statistical distributions associated with amorphous mod- 
els agree overall with the respective broadened distribu- 
tions for the counterpart crystals (see, e.g., Refs. Is^lssl) . 
which also applies to the bond-angle distributions pre- 
sented here. In addition to this, Fig.^c),(d) shows that 
the geometry optimization with DFTB results in a crys- 
tal structure that agrees with the experimental one to 
within about two percent, not only in bond distances 
and lattice constants, but also in bond angles. 

It may appear reasonable to ascribe the asymmetry of 
the distributions for S-As-S and As-S-As angles in the 
form of the shoulders on the low-angle side to the pres- 
ence of four-membered As-S-As-S rings, as has been ob- 
served in tetrahedrally bonded chalcogenide semiconduc- 
tors, e.g. GeS2^^. Although such four-membered rings 
would definitely contribute to the low-angle part of the 
distribution due to geometrical constraints, their number 
in our models is not so large (see Table ^Oj and, for that 
reason, the relative fraction of angles involved in these 
rings is rather small (about 7 % of S-As-S angles and 
13 % of As-S-As angles). We therefore conclude that the 
asymmetry of the heteropolar bond-angle distributions is 
inherent to this type of material, as is evidenced in the 
case of both crystal and amorphous structures. 

Interestingly, there are no 12-membered rings in our 
models of a-As2S3, while only such rings exist in the crys- 
tal structure of orpiment. This is another strong piece of 
evidence that models 1-4 do not contain any memory of 
the initial crystalline atomic arrangement. Also there are 
no three-membered rings, that were previously reported 
to be found in a model of a-As2Se3®, whose presence 
would contribute a few small angles and would increase 
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ring size n 
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TABLE III: Ring statistics - number of n-membered shortest 
rings for n < 22. Columns containing only zeros are 
not included. 
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FIG. 5: Fragment of model 1: two bond-sharing five- 
membered rings and the two AsSs groups connected to this 
structure. The shading of the As atoms (all with three neigh- 
bors) is darker than that of the S atoms (all with two neigh- 
bors). The HOMO level is mostly localized on the two S 
atoms that are marked black. The distance between these 
two atoms is 3.42 A. 



the number of homopolar bonds. 

Special significance can be attributed to the presence of 
five-membered rings in models 1 and 2 with an apprecia- 
ble concentration of homopolar bonds (models with all- 
heteropolar bonds contain only an even number of atoms 
in all rings). When such rings share some of the bonds, 
the resulting local structure is close to that of cage-like 
molecules (e.g. AS4S4 or AS4S3), as found in the vapor 
phase and in some chalcogenide molecular crystals. Fig[Sl 
(cf. Fig. 8.8(b) in Ref. ISo) shows two such bond-sharing 
rings. Upon breaking the two bonds connecting the rings 
to the rest of the network, the distance of 4.54 A between 
the two freed arsenic atoms could be reduced, thus pro- 
ducing another As- As homopolar bond and this group of 
atoms would then form an AS4S4 molecule. Evidence of 
the presence of such molecules in bulk AsxSi_x glasses 
from Raman-scattering experiments has recently been re- 
ported in Ref. 'st'. Our result shows that the AS4S4 frag- 
ments may not only form discrete cage- like molecules but 
also be embedded into the amorphous network. We veri- 
fied that the vibrational signatures of the AS4S4 fragment 
from models 1 and 2 are similar to those from an isolated 
AS4S4 molecule, apart from a few very symmetric modes 
of the latter. 



FIG. 6; Vibrational densities of states (a) and inverse partic- 
ipation ratios (b),(c) for models 1 and 4. The experimental 
data in (a) are obtained from Ref. l3^ 



B. Vibrational properties 

The vibrational density of states (VDOS) for models 1 
and 4 is shown in Fig. Ela). It has the two-band form 
generally observed in amorphous semiconductors. The 
vibrational spectrum for our models is essentially super- 
imposable on the calculated VDOS^ for a model of a- 
As2Se3 if the energy in Fig.El^a) is downscaled by a factor 
of about 0.67. All main features — the position of the 
gap between the acoustic and optic bands, as well as the 
relative width and height of VDOS within these bands — 
agree with available inelastic neutron-scattering experi- 
mental dataii^. Note that the experimental curve in 
Fig. EJa) corresponds to a measurement at room tem- 
perature, whereas the results obtained from simulation 
are calculated for an energy-minimum configuration in 
the harmonic approximation. An attempt to measure 
the VDOS of a-As2S3 at temperatures as low as 25 K 
was made in Ref. ITTl and the tendency for the narrowing 
and heightening of the optic band and the flattening of 
the top of the acoustic band was captured, although the 
experimental uncertainty was rather large. 

The inverse participation ratios (IPR) (see 
Fig. EIb),(c)) show that the vibrational eigenmodes 
are significantly localized in the gap between the 
acoustic and optic bands, and the high-i? end of the 
spectrum. The first three highest energy modes in 
model 1 are localized on S-S bonds, while the next 
two (with the largest IPR) are localized on IVAPs. 
The VDOS's for models 1 and 4 differ mainly in the 
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FIG. 7: Local and total electronic density of states (a) and 
inverse participation ratios (b) for model 1. The Fermi energy 
is at the orig in. The experimental data in (a) are obtained 
from Ref. IM 



absence of the just-mentioned highest energy modes 
from the spectrum corresponding to the stoichiometric 
model 4. The highest energy mode {E = 50.5 meV) in 
model 4 is localized on a relatively complex structure 
involving three AsSs pyramids in a chain As-S-As-S-As, 
where both As-S-As angles (118 and 122°) are at the 
large-0 end of the bond-angle distribution shown in 
Fig. ^d). Within the optic band, the modes with IPR 
greater than 0.2, at £^ = 43.8 meV in model 1 and at 
E = 42.9 and 45.7 meV in model 4, are localized at 
four-membered rings. In the band gap, the mode at the 
top of the acoustic band in model 1 is localized on an 
As-As-S-S chain, and the mode at the bottom of the 
optic band is localized on an IVAP. The mode at the 
top of the acoustic band in model 4 is predominantly 
localized on a four-membered ring, and the mode at the 
bottom of the optic band is almost entirely localized on 
a six-membered ring. 



C. Electronic structure 

The electronic density of states for model 1, as well 
as the inverse participation ratios for these states, are 




FIG. 8: Electronic densities of states of the models of a-As2S3 
and orpiment. 



shown in Fig. □ The calculated total EDOS IS m good 
agreement with the density of valence states2^ measured 
by X-ray photoemission spectroscopy. The total EDOS 
is also very similar to that of arsenic selenide in Ref. Hi 
where an experimental result for that material was pre- 
sented to be in agreement with the calculated one. 

The local EDOS's for different elements and orbital 
types, shown in Fig.[7{a), confirm the analysis presented 
in Ref. |^ The top of the valence band is due to the 
non-bonding lone-pair p orbitals of the S atoms, and the 
rest of the valence band is attributed to the bonding p 
orbitals on the S and As atoms. 

The s band is composed of two sub-bands. The low-i? 
sub-band at about -(15-12) eV is essentially an s type 
sulphur band, and the sub-band at about -(12-8) eV is 
predominately due to the arsenic s orbitals. In Ref. H, 
the fact that the s band of selenium is below the s band 
of arsenic was attributed to the greater nuclear charge of 
Se. Since the electronic structure of arsenic sulphide, a 
compound containing a much lighter chalcogen, virtually 
coincides with that for arsenic selenide, the above expla- 
nation is incorrect. A strong repulsion between As and 
S/Se s levels due to chemical ordering is more likely to 
be responsible for the separation of the s band into high 
(As) and low (S/Se) sub-bands^S. Note that the two s 
sub-bands are perfectly separated only when the chemi- 
cal ordering is perfect, i.e. in models with all-heteropolar 
bonds, as seen in Fig.|SId,e). The degree of admixture of 
As and S s orbitals within the low and high sub-bands, 
respectively, is about 30 % (see the "As, s" and "S, s" 
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FIG. 9: Inverse participation ratios at the band edges for all 
the four models. The vertical dashed lines in (a) demarcate 
the levels mostly localized at IVAPs, with each line type cor- 
responding to the same IVAP. 



model 


1 


2 


3 


4 


crystal 


£;(LUMO)-£(HOMO) 


1.60 


1.80 


2.00 


1.99 


2.51 


half-maximum gap 


2.67 


2.72 


2.74 


2.74 


3.13 



TABLE IV: Band-gap energies in eV estimated as the dif- 
ferences between the LUMO and HOMO eigenvalues and as 
the differences of the band-gap edges at the level 0.5 eV~^ 
in EDOS normalized so that the maximum value within the 
valence band is equal to unity. 



shown in Fig. |2| If topological (or coordination) defects 
are present in a model, as in model 1, some of the states, 
especially at the bottom of the conduction band, are lo- 
calized at them (see Fig. |^a)). This tendency was also 
emphasized in Ref. A comparable, or even greater de- 
gree of localization, is observed in model 2 (see Fig.lSJb)) 
where there are no coordination defects, but there is an 
appreciable concentration of homopolar bonds (or chem- 
ical defects). As the number of chemical defects goes to 
zero from model 2 to model 4, it is seen that localiza- 
tion at the top of the valence band remains qualitatively 
similar, whereas at the bottom of the conduction band it 
becomes less pronounced (see Fig. |5Ib)-(d)). 

The general picture is that, at the top of the va- 
lence band, the eigenstates are predominantly localized 
at what can be called sulphur-rich regions, where several 
sulphur atoms are closer than about 3.45 A, i.e. their 
interatomic distances are on the low-r side of the sec- 
ond peak in g^sii") shown in Fig. ^c) or some of these 
atoms form homopolar S-S bonds. For instance, most of 
the HOMO level in model 1 is localized at two sulphur 
atoms separated by 3.42 A and which are part of the 
molecule-like fragment depicted in Fig. [S] By inspect- 
ing the projected (local) IPRs in Fig.[7{b) at the optical 
gap edges, it is seen that the IPRs are greatest for the 
S atoms. It appears that the localization at the top of 
the valence band is facilitated by the proximity of the 
lone-pair p orbitals in the sulphur-rich regions. 

At the bottom of the conduction band, the states tend 
to localize at four-mcmbered rings in all models, and 
at S-S homopolar bonds (some of these bonds are in 
five-membered rings) and IVAPs when such defects are 
present. In model 4, all three conduction-band states 
with an IPR greater than 0.025 (see Fig. I^Jd)) are local- 
ized at four-membered rings. 



panels in Fig.[7fa)), again in agreement with Ref. l39l 

The conduction band is composed of about equal con- 
tributions from the antibonding As and S p orbitals and 
the S d orbitals. Perhaps, if the As d orbitals were in- 
cluded in the basis set, there would be a contribution 
from them too. While this might not significantly affect 
the structural and ground-state electronic properties of 
our models, this remark could be of greater importance 
for excited-state simulations. 

The EDOS's for models 1-4 and for the crystal struc- 
ture of orpiment are plotted in Fig. |H1 It is seen that 
the overall similarity is preserved for all these structures. 
The main differences among the EDOS's corresponding 
to the amorphous structures are the widening of the op- 
tical band gap (see also Table IIV|I and the clearing of 
the gap within the low-energy s-band as the number of 
defects diminishes. 

Localization of the electronic states near the optical 
band-gap edges is of great interest for studies of photoin- 
duced phenomena. A close-up of the inverse participation 
ratios at energies near the band edges for models 1-4 is 



IV. CONCLUSION 

We have generated several models of amorphous ar- 
senic sulphide by using a density-functional-based tight- 
binding method. All models agree very well with the 
neutron-scattering experimental structural data. We ob- 
serve a tendency for formation of quasi-molecular struc- 
tural groups which suggests that amorphous chalco- 
genides can be viewed as nanostructured materials. Vi- 
brational properties are also in agreement with experi- 
mental results. 

In models containing both homopolar bonds and topo- 
logical defects, a significant degree of electronic-state lo- 
calization has been observed near both band-gap edges. 
Although the coordination-number defects are optically 
active, their presence may not be necessary for exhibiting 
photostructural changes when there is a sufficient concen- 
tration of homopolar bonds in the material. This state- 
ment is supported by the observation that, upon removal 
of the coordination defects from the system, the degree of 
electronic-state localization is not reduced in the resul- 
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tant continuous network model with homopolar bonds. 
Furthermore, the valence-alternation defect concentra- 
tion is estimated to be rather small (10^^ cm~^ in Ref. fiol) 
compared to the atomic density of about 2 x 10^^ cm~'^. 
This indicates that, in the volume occupied by our 200- 
atom models, it is less likely to find a topological defect 
than to come across none. 

A stoichiometric continuous network model has al- 
lowed us to identify the structural motifs where elec- 
tronic eigenstates predominantly localize at the optical 
band edges, in the absence of coordination and chemi- 
cal defects. These are sulphur-rich regions for the top of 
the valence band and four-membered rings for the bot- 
tom of the conduction band. Electronic properties of this 
glass model with all-heteropolar bonds are very similar 
to those of the corresponding crystalline phase, orpiment, 
most notably the clean gaps in both s and p bands. Al- 
though the valence band in this case contains about as 
many localized states as in models with defects, the con- 
duction band has very few of them. It is expected that, 
if there were no four-membered rings in this structure, 
all states in the conduction band would be virtually de- 
localized. 



Therefore we conclude that perhaps the dominant con- 
tribution to photo-induced effects originates from the 
presence of electronic states localized in the vicinity of 
homopolar bonds, in support of the theoretical mod- 
els where the photo-induced structural changes are at- 
tributed to the presence of homopolar bonds in these ma- 
terials (see, e.g., Refs. I3 l4l|) . Although electronic states 
can also localize in all-heteropolar regions and in the 
vicinity of the topological defects, the contribution of 
such states is likely to be rather small due to the low 
degree of localization in the conduction band and the 
low concentration of such defects, respectively. Verifica- 
tion of this conjecture requires excited-state calculations 
and is beyond the scope of the present paper. We plan 
to do these calculations in the future. 
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